R Notebook
package 㤼㸱tidyverse㤼㸲 was built under R version 3.6.3package 㤼㸱ggplot2㤼㸲 was built under R version 3.6.3package 㤼㸱tibble㤼㸲 was built under R version 3.6.3package 㤼㸱tidyr㤼㸲 was built under R version 3.6.3package 㤼㸱purrr㤼㸲 was built under R version 3.6.3package 㤼㸱dplyr㤼㸲 was built under R version 3.6.3package 㤼㸱ReacTran㤼㸲 was built under R version 3.6.3package 㤼㸱rootSolve㤼㸲 was built under R version 3.6.3package 㤼㸱deSolve㤼㸲 was built under R version 3.6.3package 㤼㸱shape㤼㸲 was built under R version 3.6.3Woah! Welcome back to having to think after the collective food coma of the Thanksgiving break! Remember the good times we had talking about how the atmospheric signal of water vapor got incoded into the ice sheet? If you don’t, here is a quick reminder:
A little reminder of where we are in the system.
We are in this “viscous sublayer” or VSL (also called the laminar layer in the Craig-Gordon model). In this layer, transport of different isotopes is entirely diffusion! In these exercises we are going to explore what the means for the top of the ice sheet (or “firn”). More concisely, we are going to use a toy diffusion model to see the effect that the VSL’s thickness has on what the firn “sees” from the atmosphere. From the paper:
“A positive linear relation between B_best and delta_z is expected since the influence of the bottom boundary condition is more attenuated for a thicker VSL.”
What the hell is B_best anyways? Well we have snow measurements and atmospheric measurements but it’s hard to make a measurement AT the snow surface. Madsen et al. assume that there is some dirunal variation B, a mean value A, and phase offset C. In the end, the value at the snow surface will look something like this:
\[ \delta^*(z=0,t) = A^* + B^*sin(\frac{2\pi t}{T_{day}}+ C^*) \] Exercise 1: Go ahead and make a function called “dirurnal_isotope” as a function of A, B, C, and t. Then lets do a sanity check by plotting it in a tibble called “isotope_model_data”. Go ahead and use the values from the 26th-29th of June from the paper for \[\delta^{18}O\] which so happens to be \[A_{best} = -43.9, B_{best} = 3.3, C_{best} = 14.4\].
diurnal_isotope <- function(A, B, C, t) {
isotope <- A + B*sin((2*pi*t)/24 + C)
return(isotope)
}Here is the plot.
isotope_model_data <- tibble(time = seq(from = 0, to = 24, by = 1), snow_value = diurnal_isotope(-43.9, 3.3, 14.4, time), air_value = diurnal_isotope(-45, 10, 14.0, time))
isotope_model_data %>%
ggplot(aes(x = time, y = snow_value)) +
geom_line(lwd = 1) +
theme(text = element_text(size = 20))
Differential equationsin R isn’t crazy mature, though it has some great basic packages that will help us get pretty close. We are going to solve the diffusion equation in 1D with the help of ReacTran (which will give us an equation) and deSolve (which will do our solving). We will get a steady state solution at each time step, go much more forward in time and solve again with the new atmoshperic boundary conditions. The very first thing to do is get comfortable with solving the diffusion equation and inform our model spin-up time.
Note: We are ignoring the advection term here. Though Madesen et al. includes it, it is pretty damn complicated and outside the scope of these excersies.
# Lets make a grid in the spatial dimensions. This will represent a VSL of 0.1 m
Grid <- setup.grid.1D(N = 10000, L = 0.01)
# The function for solving
pde1D <-function(t, C, parms) {
tran <- tran.1D(C = C, D = D,
C.down = Cdown, C.up = Cup, dx = Grid)$dC
list(tran) # return value: rate of change
}
D <- 0.00000009723 # diffusion constant for 18O in 0.1m/sec
Cdown <- diurnal_isotope(-43.9, 3.3, 14.4, 0) #snow value
Cup <- diurnal_isotope(-45, 10, 14.4, 0) #air value
Cave <- (Cdown + Cup) / 2 # we need an inital condition somehow, lets do the average
diff_time <- 180 #number of time steps in units of what our diffusion constant is in time
times <- seq(0, diff_time, by = 1) #make time array
system.time(
out <- ode.1D(y = rep(Cave, Grid$N),
times = times, func = pde1D,
parms = NULL, nspec = 1)
) user system elapsed
1.47 0.02 1.61 tail(out[, 1:10], n = 10) #sample the out file from the solver time 1 2 3 4 5 6 7 8 9
[172,] 171 -35.34369 -35.34423 -35.34477 -35.34531 -35.34585 -35.34638 -35.34692 -35.34746 -35.348
[173,] 172 -35.34369 -35.34423 -35.34477 -35.34531 -35.34585 -35.34638 -35.34692 -35.34746 -35.348
[174,] 173 -35.34369 -35.34423 -35.34477 -35.34531 -35.34585 -35.34638 -35.34692 -35.34746 -35.348
[175,] 174 -35.34369 -35.34423 -35.34477 -35.34531 -35.34584 -35.34638 -35.34692 -35.34746 -35.348
[176,] 175 -35.34369 -35.34423 -35.34477 -35.34531 -35.34584 -35.34638 -35.34692 -35.34746 -35.348
[177,] 176 -35.34369 -35.34423 -35.34477 -35.34531 -35.34584 -35.34638 -35.34692 -35.34746 -35.348
[178,] 177 -35.34369 -35.34423 -35.34477 -35.34531 -35.34584 -35.34638 -35.34692 -35.34746 -35.348
[179,] 178 -35.34369 -35.34423 -35.34477 -35.34531 -35.34584 -35.34638 -35.34692 -35.34746 -35.348
[180,] 179 -35.34369 -35.34423 -35.34477 -35.34531 -35.34584 -35.34638 -35.34692 -35.34746 -35.348
[181,] 180 -35.34369 -35.34423 -35.34477 -35.34531 -35.34584 -35.34638 -35.34692 -35.34746 -35.348Let’s take the solution and plot it. Also we’ll make sure to take the last time file as our future \[t = 0\] boundary condition and call it steady_state_solution.
sss_temp <- out[121,]
steady_state_solution <- sss_temp[-1]
image(out, xlab = "time, sec",
ylab = "Distance, m",
main = "delta 18O at the start", add.contour = TRUE)
To do this at each time step we’ll pack this into one big function, “isotope_vapor_diff”. We just care about the end values so that’s the output. For the data book-keeping, tibbles can’t hold multi-dimenional data at a point unless it is a list. I have packed away the output as a list and unpacked it from the “starting_values” input variable.
isotope_vapor_diff = function(starting_values, snow_value, atmos_value) {
if(is.list(starting_values)){
starting_values <- unlist(starting_values)
}
D <- 0.00000009723 # diffusion constant
Cdown <- snow_value
Cup <- atmos_value
diff_time <- 600
Grid <- setup.grid.1D(N = 10000, L = 0.01)
pde1D <-function(t, C, parms) {
tran <- tran.1D(C = C, D = D,
C.down = Cdown, C.up = Cup, dx = Grid)$dC
list(tran) # return value: rate of change
}
times <- seq(0, diff_time, by = 1)
out <- ode.1D(y = starting_values,
times = times, func = pde1D,
parms = NULL, nspec = 1)
sss_temp <- out[121,]
steady_state_solution <- list(sss_temp[-1])
return(steady_state_solution)
}Exercise 2: Use this fancy new function to explore how the VSL would look in three minutes if you changed the surface snow value and air value from our previous steady state solution to two new values. Plot your results.
Note: Remember that the function returns a list. You can unpack that bad boy with a statement that looks like “unlist(your_list_here)”.
func_test <- isotope_vapor_diff(steady_state_solution, -45, -40)
plot(seq(1,10000), unlist(func_test))
On to the main event. Now we need to add in our first steady state solution to serve as the t = 0 boundary of our model!
isotope_model_data <- isotope_model_data %>%
add_column(isotope_space = NA)
isotope_model_data$isotope_space[1] = list(steady_state_solution)
isotope_model_dataLooks good. Now we use the diffusion function to fill out the rest of the isotope_space column using the end of previous step as the starting value.
interator = length(isotope_model_data$isotope_space) - 1
system.time(
for (i in 1:interator) {
isotope_model_data$isotope_space[i+1] = isotope_vapor_diff(isotope_model_data$isotope_space[i], isotope_model_data$snow_value[i], isotope_model_data$air_value[i])
}
) user system elapsed
45.16 0.28 45.55 isotope_model_dataThis is how we might look at a single solution. In this case, hour 15.
sample_from_isotope_tibble <- unlist(isotope_model_data$isotope_space[4])
plot(seq(1,10000), sample_from_isotope_tibble)
At each time, we can look at the isotope value at different parts of the space as a rough approximation of what value the snow would see if the VSL was that thick.
Exercise 3: Add to the tibble four new columns for the value at four different heights of your choosing. Is your result linear with respect to VSL height? How could you tell?
isotope_model_data <- isotope_model_data %>% rowwise() %>%
mutate(value_.2m = unlist(isotope_space)[8000])
isotope_model_data <- isotope_model_data %>% rowwise() %>%
mutate(value_.4m = unlist(isotope_space)[6000])
isotope_model_data <- isotope_model_data %>% rowwise() %>%
mutate(value_.6m = unlist(isotope_space)[4000])
isotope_model_data <- isotope_model_data %>% rowwise() %>%
mutate(value_.8m = unlist(isotope_space)[2000])
isotope_model_datadiff_heights <- isotope_model_data %>%
ggplot(aes(x = time)) +
geom_line(aes(y = value_.2m), lwd = 2) +
geom_line(aes(y = value_.4m), color="blue", lwd = 2) +
geom_line(aes(y = value_.6m), color="green", lwd = 2) +
geom_line(aes(y = value_.8m), color="cyan", lwd = 2) +
theme(text = element_text(size = 20))
diff_heights
LS0tDQp0aXRsZTogIlIgTm90ZWJvb2siDQoNCm91dHB1dDoNCiAgcGRmX2RvY3VtZW50OiBkZWZhdWx0DQogIGZpZ193aWR0aDogMTIgDQogIGZpZ19oZWlnaHQ6IDggDQotLS0NCmBgYHtyLCBtZXNzYWdlPUZBTFNFLCBlY2hvPUZBTFNFfQ0KbGlicmFyeSh0aWR5dmVyc2UpDQpsaWJyYXJ5KFJlYWNUcmFuKQ0KcmVxdWlyZShkZVNvbHZlKQ0KYGBgDQoNCldvYWghIFdlbGNvbWUgYmFjayB0byBoYXZpbmcgdG8gdGhpbmsgYWZ0ZXIgdGhlIGNvbGxlY3RpdmUgZm9vZCBjb21hIG9mIHRoZSBUaGFua3NnaXZpbmcgYnJlYWshIFJlbWVtYmVyIHRoZSBnb29kIHRpbWVzIHdlIGhhZCB0YWxraW5nIGFib3V0IGhvdyB0aGUgYXRtb3NwaGVyaWMgc2lnbmFsIG9mIHdhdGVyIHZhcG9yIGdvdCBpbmNvZGVkIGludG8gdGhlIGljZSBzaGVldD8gSWYgeW91IGRvbid0LCBoZXJlIGlzIGEgcXVpY2sgcmVtaW5kZXI6DQoNCiFbQSBsaXR0bGUgcmVtaW5kZXIgb2Ygd2hlcmUgd2UgYXJlIGluIHRoZSBzeXN0ZW0uXShtYWRzZW5fdG93ZXIuUE5HKQ0KDQpXZSBhcmUgaW4gdGhpcyAidmlzY291cyBzdWJsYXllciIgb3IgVlNMIChhbHNvIGNhbGxlZCB0aGUgbGFtaW5hciBsYXllciBpbiB0aGUgQ3JhaWctR29yZG9uIG1vZGVsKS4gSW4gdGhpcyBsYXllciwgdHJhbnNwb3J0IG9mIGRpZmZlcmVudCBpc290b3BlcyBpcyBlbnRpcmVseSBkaWZmdXNpb24hIEluIHRoZXNlIGV4ZXJjaXNlcyB3ZSBhcmUgZ29pbmcgdG8gZXhwbG9yZSB3aGF0IHRoZSBtZWFucyBmb3IgdGhlIHRvcCBvZiB0aGUgaWNlIHNoZWV0IChvciAiZmlybiIpLiBNb3JlIGNvbmNpc2VseSwgd2UgYXJlIGdvaW5nIHRvIHVzZSBhIHRveSBkaWZmdXNpb24gbW9kZWwgdG8gc2VlIHRoZSBlZmZlY3QgdGhhdCB0aGUgVlNMJ3MgdGhpY2tuZXNzIGhhcyBvbiB3aGF0IHRoZSBmaXJuICJzZWVzIiBmcm9tIHRoZSBhdG1vc3BoZXJlLiBGcm9tIHRoZSBwYXBlcjoNCg0KKuKAnEEgcG9zaXRpdmUgbGluZWFyIHJlbGF0aW9uIGJldHdlZW4gQl9iZXN0IGFuZCBkZWx0YV96IGlzIGV4cGVjdGVkIHNpbmNlIHRoZSBpbmZsdWVuY2Ugb2YgdGhlIGJvdHRvbSBib3VuZGFyeSBjb25kaXRpb24gaXMgbW9yZSBhdHRlbnVhdGVkIGZvciBhIHRoaWNrZXIgVlNMLuKAnSoNCg0KV2hhdCB0aGUgaGVsbCBpcyBCX2Jlc3QgYW55d2F5cz8gV2VsbCB3ZSBoYXZlIHNub3cgbWVhc3VyZW1lbnRzIGFuZCBhdG1vc3BoZXJpYyBtZWFzdXJlbWVudHMgYnV0IGl0J3MgaGFyZCB0byBtYWtlIGEgbWVhc3VyZW1lbnQgKipBVCoqIHRoZSBzbm93IHN1cmZhY2UuIE1hZHNlbiBldCBhbC4gYXNzdW1lIHRoYXQgdGhlcmUgaXMgc29tZSBkaXJ1bmFsIHZhcmlhdGlvbiBCLCBhIG1lYW4gdmFsdWUgQSwgYW5kIHBoYXNlIG9mZnNldCBDLiBJbiB0aGUgZW5kLCB0aGUgdmFsdWUgYXQgdGhlIHNub3cgc3VyZmFjZSB3aWxsIGxvb2sgc29tZXRoaW5nIGxpa2UgdGhpczoNCg0KJCQNClxkZWx0YV4qKHo9MCx0KSA9IEFeKiArIEJeKnNpbihcZnJhY3syXHBpIHR9e1Rfe2RheX19KyBDXiopDQokJA0KKipFeGVyY2lzZSAxOioqIEdvIGFoZWFkIGFuZCBtYWtlIGEgZnVuY3Rpb24gY2FsbGVkICJkaXJ1cm5hbF9pc290b3BlIiBhcyBhIGZ1bmN0aW9uIG9mIEEsIEIsIEMsIGFuZCB0LiBUaGVuIGxldHMgZG8gYSBzYW5pdHkgY2hlY2sgYnkgcGxvdHRpbmcgaXQgaW4gYSB0aWJibGUgY2FsbGVkICJpc290b3BlX21vZGVsX2RhdGEiLiBHbyBhaGVhZCBhbmQgdXNlIHRoZSB2YWx1ZXMgZnJvbSB0aGUgMjZ0aC0yOXRoIG9mIEp1bmUgZnJvbSB0aGUgcGFwZXIgZm9yICQkXGRlbHRhXnsxOH1PJCQgd2hpY2ggc28gaGFwcGVucyB0byBiZSAkJEFfe2Jlc3R9ID0gLTQzLjksIEJfe2Jlc3R9ID0gMy4zLCBDX3tiZXN0fSA9IDE0LjQkJC4gDQoNCmBgYHtyfQ0KZGl1cm5hbF9pc290b3BlIDwtIGZ1bmN0aW9uKEEsIEIsIEMsIHQpIHsNCiAgIGlzb3RvcGUgPC0gQSArIEIqc2luKCgyKnBpKnQpLzI0ICsgQykNCiAgIHJldHVybihpc290b3BlKQ0KfQ0KYGBgDQoNCkhlcmUgaXMgdGhlIHBsb3QuDQoNCmBgYHtyLCBmaWcud2lkdGg9MTAsIGZpZy5oZWlnaHQ9NSwgZHBpPTYwMH0NCmlzb3RvcGVfbW9kZWxfZGF0YSA8LSB0aWJibGUodGltZSA9IHNlcShmcm9tID0gMCwgdG8gPSAyNCwgYnkgPSAxKSwgc25vd192YWx1ZSA9IGRpdXJuYWxfaXNvdG9wZSgtNDMuOSwgMy4zLCAxNC40LCB0aW1lKSwgYWlyX3ZhbHVlID0gZGl1cm5hbF9pc290b3BlKC00NSwgMTAsIDE0LjAsIHRpbWUpKQ0KDQppc290b3BlX21vZGVsX2RhdGEgJT4lDQogICAgZ2dwbG90KGFlcyh4ID0gdGltZSwgeSA9IHNub3dfdmFsdWUpKSArDQogICAgZ2VvbV9saW5lKGx3ZCA9IDEpICsgDQogICAgdGhlbWUodGV4dCA9IGVsZW1lbnRfdGV4dChzaXplID0gMjApKQ0KYGBgDQpEaWZmZXJlbnRpYWwgZXF1YXRpb25zaW4gUiBpc24ndCBjcmF6eSBtYXR1cmUsIHRob3VnaCBpdCBoYXMgc29tZSBncmVhdCBiYXNpYyBwYWNrYWdlcyB0aGF0IHdpbGwgaGVscCB1cyBnZXQgcHJldHR5IGNsb3NlLiBXZSBhcmUgZ29pbmcgdG8gc29sdmUgdGhlIGRpZmZ1c2lvbiBlcXVhdGlvbiBpbiAxRCB3aXRoIHRoZSBoZWxwIG9mIFJlYWNUcmFuICh3aGljaCB3aWxsIGdpdmUgdXMgYW4gZXF1YXRpb24pIGFuZCBkZVNvbHZlICh3aGljaCB3aWxsIGRvIG91ciBzb2x2aW5nKS4gV2Ugd2lsbCBnZXQgYSBzdGVhZHkgc3RhdGUgc29sdXRpb24gYXQgZWFjaCB0aW1lIHN0ZXAsIGdvIG11Y2ggbW9yZSBmb3J3YXJkIGluIHRpbWUgYW5kIHNvbHZlIGFnYWluIHdpdGggdGhlIG5ldyBhdG1vc2hwZXJpYyBib3VuZGFyeSBjb25kaXRpb25zLiBUaGUgdmVyeSBmaXJzdCB0aGluZyB0byBkbyBpcyBnZXQgY29tZm9ydGFibGUgd2l0aCBzb2x2aW5nIHRoZSBkaWZmdXNpb24gZXF1YXRpb24gYW5kIGluZm9ybSBvdXIgbW9kZWwgc3Bpbi11cCB0aW1lLg0KDQoqKk5vdGU6KiogV2UgYXJlIGlnbm9yaW5nIHRoZSBhZHZlY3Rpb24gdGVybSBoZXJlLiBUaG91Z2ggTWFkZXNlbiBldCBhbC4gaW5jbHVkZXMgaXQsIGl0IGlzIHByZXR0eSBkYW1uIGNvbXBsaWNhdGVkIGFuZCBvdXRzaWRlIHRoZSBzY29wZSBvZiB0aGVzZSBleGNlcnNpZXMuDQoNCmBgYHtyfQ0KIyBMZXRzIG1ha2UgYSBncmlkIGluIHRoZSBzcGF0aWFsIGRpbWVuc2lvbnMuIFRoaXMgd2lsbCByZXByZXNlbnQgYSBWU0wgb2YgMC4xIG0NCkdyaWQgPC0gc2V0dXAuZ3JpZC4xRChOID0gMTAwMDAsIEwgPSAwLjAxKQ0KDQojIFRoZSBmdW5jdGlvbiBmb3Igc29sdmluZw0KcGRlMUQgPC1mdW5jdGlvbih0LCBDLCBwYXJtcykgew0KICB0cmFuIDwtIHRyYW4uMUQoQyA9IEMsIEQgPSBELA0KICBDLmRvd24gPSBDZG93biwgQy51cCA9IEN1cCwgZHggPSBHcmlkKSRkQw0KICBsaXN0KHRyYW4pICMgcmV0dXJuIHZhbHVlOiByYXRlIG9mIGNoYW5nZQ0KfQ0KDQpEIDwtIDAuMDAwMDAwMDk3MjMgIyBkaWZmdXNpb24gY29uc3RhbnQgZm9yIDE4TyBpbiAwLjFtL3NlYw0KQ2Rvd24gPC0gZGl1cm5hbF9pc290b3BlKC00My45LCAzLjMsIDE0LjQsIDApICNzbm93IHZhbHVlDQpDdXAgPC0gZGl1cm5hbF9pc290b3BlKC00NSwgMTAsIDE0LjQsIDApICNhaXIgdmFsdWUNCkNhdmUgPC0gKENkb3duICsgQ3VwKSAvIDIgIyB3ZSBuZWVkIGFuIGluaXRhbCBjb25kaXRpb24gc29tZWhvdywgbGV0cyBkbyB0aGUgYXZlcmFnZQ0KZGlmZl90aW1lIDwtIDE4MCAjbnVtYmVyIG9mIHRpbWUgc3RlcHMgaW4gdW5pdHMgb2Ygd2hhdCBvdXIgZGlmZnVzaW9uIGNvbnN0YW50IGlzIGluIHRpbWUNCg0KdGltZXMgPC0gc2VxKDAsIGRpZmZfdGltZSwgYnkgPSAxKSAjbWFrZSB0aW1lIGFycmF5DQpzeXN0ZW0udGltZSgNCiAgb3V0IDwtIG9kZS4xRCh5ID0gcmVwKENhdmUsIEdyaWQkTiksDQogIHRpbWVzID0gdGltZXMsIGZ1bmMgPSBwZGUxRCwNCiAgcGFybXMgPSBOVUxMLCBuc3BlYyA9IDEpDQopDQoNCnRhaWwob3V0WywgMToxMF0sIG4gPSAxMCkgI3NhbXBsZSB0aGUgb3V0IGZpbGUgZnJvbSB0aGUgc29sdmVyDQoNCmBgYA0KDQpMZXQncyB0YWtlIHRoZSBzb2x1dGlvbiBhbmQgcGxvdCBpdC4gQWxzbyB3ZSdsbCBtYWtlIHN1cmUgdG8gdGFrZSB0aGUgbGFzdCB0aW1lIGZpbGUgYXMgb3VyIGZ1dHVyZSAkJHQgPSAwJCQgYm91bmRhcnkgY29uZGl0aW9uIGFuZCBjYWxsIGl0IHN0ZWFkeV9zdGF0ZV9zb2x1dGlvbi4gDQoNCmBgYHtyfQ0Kc3NzX3RlbXAgPC0gb3V0WzEyMSxdDQpzdGVhZHlfc3RhdGVfc29sdXRpb24gPC0gc3NzX3RlbXBbLTFdDQppbWFnZShvdXQsIHhsYWIgPSAidGltZSwgc2VjIiwNCiAgeWxhYiA9ICJEaXN0YW5jZSwgbSIsDQogIG1haW4gPSAiZGVsdGEgMThPIGF0IHRoZSBzdGFydCIsIGFkZC5jb250b3VyID0gVFJVRSkNCmBgYA0KDQpUbyBkbyB0aGlzIGF0IGVhY2ggdGltZSBzdGVwIHdlJ2xsIHBhY2sgdGhpcyBpbnRvIG9uZSBiaWcgZnVuY3Rpb24sICJpc290b3BlX3ZhcG9yX2RpZmYiLiBXZSBqdXN0IGNhcmUgYWJvdXQgdGhlIGVuZCB2YWx1ZXMgc28gdGhhdCdzIHRoZSBvdXRwdXQuIEZvciB0aGUgZGF0YSBib29rLWtlZXBpbmcsIHRpYmJsZXMgY2FuJ3QgaG9sZCBtdWx0aS1kaW1lbmlvbmFsIGRhdGEgYXQgYSBwb2ludCB1bmxlc3MgaXQgaXMgYSBsaXN0LiBJIGhhdmUgcGFja2VkIGF3YXkgdGhlIG91dHB1dCBhcyBhIGxpc3QgYW5kIHVucGFja2VkIGl0IGZyb20gdGhlICJzdGFydGluZ192YWx1ZXMiIGlucHV0IHZhcmlhYmxlLg0KDQpgYGB7cn0NCmlzb3RvcGVfdmFwb3JfZGlmZiA9IGZ1bmN0aW9uKHN0YXJ0aW5nX3ZhbHVlcywgc25vd192YWx1ZSwgYXRtb3NfdmFsdWUpIHsNCiAgaWYoaXMubGlzdChzdGFydGluZ192YWx1ZXMpKXsNCiAgICBzdGFydGluZ192YWx1ZXMgPC0gdW5saXN0KHN0YXJ0aW5nX3ZhbHVlcykNCiAgfQ0KICANCiAgRCA8LSAwLjAwMDAwMDA5NzIzICMgZGlmZnVzaW9uIGNvbnN0YW50DQogIENkb3duIDwtIHNub3dfdmFsdWUNCiAgQ3VwIDwtIGF0bW9zX3ZhbHVlDQogIGRpZmZfdGltZSA8LSA2MDANCg0KICANCiAgR3JpZCA8LSBzZXR1cC5ncmlkLjFEKE4gPSAxMDAwMCwgTCA9IDAuMDEpDQoNCiAgcGRlMUQgPC1mdW5jdGlvbih0LCBDLCBwYXJtcykgew0KICAgIHRyYW4gPC0gdHJhbi4xRChDID0gQywgRCA9IEQsDQogICAgQy5kb3duID0gQ2Rvd24sIEMudXAgPSBDdXAsIGR4ID0gR3JpZCkkZEMNCiAgICBsaXN0KHRyYW4pICMgcmV0dXJuIHZhbHVlOiByYXRlIG9mIGNoYW5nZQ0KICAgIH0NCg0KDQogIHRpbWVzIDwtIHNlcSgwLCBkaWZmX3RpbWUsIGJ5ID0gMSkNCiAgb3V0IDwtIG9kZS4xRCh5ID0gc3RhcnRpbmdfdmFsdWVzLA0KICAgIHRpbWVzID0gdGltZXMsIGZ1bmMgPSBwZGUxRCwNCiAgICBwYXJtcyA9IE5VTEwsIG5zcGVjID0gMSkNCg0KICBzc3NfdGVtcCA8LSBvdXRbMTIxLF0NCiAgc3RlYWR5X3N0YXRlX3NvbHV0aW9uIDwtIGxpc3Qoc3NzX3RlbXBbLTFdKQ0KICByZXR1cm4oc3RlYWR5X3N0YXRlX3NvbHV0aW9uKQ0KfQ0KYGBgDQoNCioqRXhlcmNpc2UgMjoqKiBVc2UgdGhpcyBmYW5jeSBuZXcgZnVuY3Rpb24gdG8gZXhwbG9yZSBob3cgdGhlIFZTTCB3b3VsZCBsb29rIGluIHRocmVlIG1pbnV0ZXMgaWYgeW91IGNoYW5nZWQgdGhlIHN1cmZhY2Ugc25vdyB2YWx1ZSBhbmQgYWlyIHZhbHVlIGZyb20gb3VyIHByZXZpb3VzIHN0ZWFkeSBzdGF0ZSBzb2x1dGlvbiB0byB0d28gbmV3IHZhbHVlcy4gUGxvdCB5b3VyIHJlc3VsdHMuDQoNCioqTm90ZToqKiBSZW1lbWJlciB0aGF0IHRoZSBmdW5jdGlvbiByZXR1cm5zIGEgbGlzdC4gWW91IGNhbiB1bnBhY2sgdGhhdCBiYWQgYm95IHdpdGggYSBzdGF0ZW1lbnQgdGhhdCBsb29rcyBsaWtlICJ1bmxpc3QoeW91cl9saXN0X2hlcmUpIi4NCg0KYGBge3J9DQpmdW5jX3Rlc3QgPC0gaXNvdG9wZV92YXBvcl9kaWZmKHN0ZWFkeV9zdGF0ZV9zb2x1dGlvbiwgLTQ1LCAtNDApDQpwbG90KHNlcSgxLDEwMDAwKSwgdW5saXN0KGZ1bmNfdGVzdCkpDQpgYGANCg0KT24gdG8gdGhlIG1haW4gZXZlbnQuIE5vdyB3ZSBuZWVkIHRvIGFkZCBpbiBvdXIgZmlyc3Qgc3RlYWR5IHN0YXRlIHNvbHV0aW9uIHRvIHNlcnZlIGFzIHRoZSB0ID0gMCBib3VuZGFyeSBvZiBvdXIgbW9kZWwhDQoNCmBgYHtyfQ0KaXNvdG9wZV9tb2RlbF9kYXRhIDwtIGlzb3RvcGVfbW9kZWxfZGF0YSAlPiUNCiAgYWRkX2NvbHVtbihpc290b3BlX3NwYWNlID0gTkEpDQppc290b3BlX21vZGVsX2RhdGEkaXNvdG9wZV9zcGFjZVsxXSA9IGxpc3Qoc3RlYWR5X3N0YXRlX3NvbHV0aW9uKQ0KaXNvdG9wZV9tb2RlbF9kYXRhDQpgYGANCg0KTG9va3MgZ29vZC4gTm93IHdlIHVzZSB0aGUgZGlmZnVzaW9uIGZ1bmN0aW9uIHRvIGZpbGwgb3V0IHRoZSByZXN0IG9mIHRoZSBpc290b3BlX3NwYWNlIGNvbHVtbiB1c2luZyB0aGUgZW5kIG9mIHByZXZpb3VzIHN0ZXAgYXMgdGhlIHN0YXJ0aW5nIHZhbHVlLg0KDQpgYGB7cn0NCmludGVyYXRvciA9IGxlbmd0aChpc290b3BlX21vZGVsX2RhdGEkaXNvdG9wZV9zcGFjZSkgLSAxDQpzeXN0ZW0udGltZSgNCmZvciAoaSBpbiAxOmludGVyYXRvcikgew0KICBpc290b3BlX21vZGVsX2RhdGEkaXNvdG9wZV9zcGFjZVtpKzFdID0gaXNvdG9wZV92YXBvcl9kaWZmKGlzb3RvcGVfbW9kZWxfZGF0YSRpc290b3BlX3NwYWNlW2ldLCBpc290b3BlX21vZGVsX2RhdGEkc25vd192YWx1ZVtpXSwgaXNvdG9wZV9tb2RlbF9kYXRhJGFpcl92YWx1ZVtpXSkNCn0NCikNCmlzb3RvcGVfbW9kZWxfZGF0YQ0KYGBgDQoNClRoaXMgaXMgaG93IHdlIG1pZ2h0IGxvb2sgYXQgYSBzaW5nbGUgc29sdXRpb24uIEluIHRoaXMgY2FzZSwgaG91ciAxNS4NCg0KYGBge3J9DQpzYW1wbGVfZnJvbV9pc290b3BlX3RpYmJsZSA8LSB1bmxpc3QoaXNvdG9wZV9tb2RlbF9kYXRhJGlzb3RvcGVfc3BhY2VbNF0pDQpwbG90KHNlcSgxLDEwMDAwKSwgc2FtcGxlX2Zyb21faXNvdG9wZV90aWJibGUpDQpgYGANCg0KQXQgZWFjaCB0aW1lLCB3ZSBjYW4gbG9vayBhdCB0aGUgaXNvdG9wZSB2YWx1ZSBhdCBkaWZmZXJlbnQgcGFydHMgb2YgdGhlIHNwYWNlIGFzIGEgcm91Z2ggYXBwcm94aW1hdGlvbiBvZiB3aGF0IHZhbHVlIHRoZSBzbm93IHdvdWxkIHNlZSBpZiB0aGUgVlNMIHdhcyB0aGF0IHRoaWNrLg0KDQoqKkV4ZXJjaXNlIDM6KiogQWRkIHRvIHRoZSB0aWJibGUgZm91ciBuZXcgY29sdW1ucyBmb3IgdGhlIHZhbHVlIGF0IGZvdXIgZGlmZmVyZW50IGhlaWdodHMgb2YgeW91ciBjaG9vc2luZy4gSXMgeW91ciByZXN1bHQgbGluZWFyIHdpdGggcmVzcGVjdCB0byBWU0wgaGVpZ2h0PyBIb3cgY291bGQgeW91IHRlbGw/DQoNCmBgYHtyfQ0KaXNvdG9wZV9tb2RlbF9kYXRhIDwtIGlzb3RvcGVfbW9kZWxfZGF0YSAlPiUgcm93d2lzZSgpICU+JQ0KICBtdXRhdGUodmFsdWVfLjJtID0gdW5saXN0KGlzb3RvcGVfc3BhY2UpWzgwMDBdKQ0KDQppc290b3BlX21vZGVsX2RhdGEgPC0gaXNvdG9wZV9tb2RlbF9kYXRhICU+JSByb3d3aXNlKCkgJT4lDQogIG11dGF0ZSh2YWx1ZV8uNG0gPSB1bmxpc3QoaXNvdG9wZV9zcGFjZSlbNjAwMF0pDQoNCmlzb3RvcGVfbW9kZWxfZGF0YSA8LSBpc290b3BlX21vZGVsX2RhdGEgJT4lIHJvd3dpc2UoKSAlPiUNCiAgbXV0YXRlKHZhbHVlXy42bSA9IHVubGlzdChpc290b3BlX3NwYWNlKVs0MDAwXSkNCg0KaXNvdG9wZV9tb2RlbF9kYXRhIDwtIGlzb3RvcGVfbW9kZWxfZGF0YSAlPiUgcm93d2lzZSgpICU+JQ0KICBtdXRhdGUodmFsdWVfLjhtID0gdW5saXN0KGlzb3RvcGVfc3BhY2UpWzIwMDBdKQ0KDQppc290b3BlX21vZGVsX2RhdGENCmBgYA0KDQpgYGB7ciwgZmlnLndpZHRoPTEwLCBmaWcuaGVpZ2h0PTUsIGRwaSA9IDMwMH0NCmRpZmZfaGVpZ2h0cyA8LSBpc290b3BlX21vZGVsX2RhdGEgJT4lDQogIGdncGxvdChhZXMoeCA9IHRpbWUpKSArIA0KICBnZW9tX2xpbmUoYWVzKHkgPSB2YWx1ZV8uMm0pLCBsd2QgPSAyKSArIA0KICBnZW9tX2xpbmUoYWVzKHkgPSB2YWx1ZV8uNG0pLCBjb2xvcj0iYmx1ZSIsIGx3ZCA9IDIpICsNCiAgZ2VvbV9saW5lKGFlcyh5ID0gdmFsdWVfLjZtKSwgY29sb3I9ImdyZWVuIiwgbHdkID0gMikgKyANCiAgZ2VvbV9saW5lKGFlcyh5ID0gdmFsdWVfLjhtKSwgY29sb3I9ImN5YW4iLCBsd2QgPSAyKSArIA0KICB0aGVtZSh0ZXh0ID0gZWxlbWVudF90ZXh0KHNpemUgPSAyMCkpDQpkaWZmX2hlaWdodHMNCmBgYA0KDQo=